Converting Ferroptosis heatmap code to notebook

Code copied from Ferroptosis_heatmaps.R and fixed to run within this repo.

library(reshape2)
library(ggplot2)
library(scales)
library(pheatmap)

SAVEPLOTS <- FALSE

load("../data/merged_countdata.RData") # object countdata
Ferr <- read.delim("../data/KEGGFerroptosis_hsa04216_06-25-18.txt", header=T, stringsAsFactors = F)
ferr_genes <- Ferr$GeneName  # ferroptosis genes from KEGG (2018)
load("../data/RLD_SC-1,7,10_0,3,8d_20180701.RData")  # object rld

Loading required package: DESeq2
Loading required package: S4Vectors
Warning: package ‘S4Vectors’ was built under R version 4.3.2
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position,
    rank, rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, union, unique, unsplit, which.max, which.min


Attaching package: ‘S4Vectors’

The following object is masked from ‘package:utils’:

    findMatches

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Loading required package: IRanges
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Warning: package ‘GenomeInfoDb’ was built under R version 4.3.2
Loading required package: SummarizedExperiment
Warning: package ‘SummarizedExperiment’ was built under R version 4.3.2
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs,
    colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds,
    colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse,
    rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads,
    rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2,
    rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for
    packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians

Another ferroptosis gene list from Wikipathways/GSEA: https://www.gsea-msigdb.org/gsea/msigdb/cards/WP_FERROPTOSIS.html File downloaded 2024-01-14 and saved in ../data/WP_FERROPTOSIS.v2023.2.Hs.tsv

ferr_genes_wikipw_raw <- read.table('../data/WP_FERROPTOSIS.v2023.2.Hs.tsv', sep="\t", row.names=1)
ferr_genes_wikipw <- ferr_genes_wikipw_raw["GENE_SYMBOLS",]
ferr_genes_wikipw <- unlist(strsplit(ferr_genes_wikipw,","))
ferr_genes_wikipw <- ferr_genes_wikipw[ferr_genes_wikipw!=""]

Unclear why these genes should be excluded.

excl_genes <- c("ALOX15", "ACSL1", "ACSL3",
                "ACSL5", "ACSL6", "TP53", "TF",
                "CP", "MAP1LC3A", "MAP1LC3C",
                "CYBB")

Define which gene set to use for all plots

object is ferrgenes_4plots. Setting to first item in list will use newer and larger list of genes from wikipathways; setting to item 2 will make output match previously generated plots.

ferr_genes_4plots <- list(ferr_genes_wikipw,ferr_genes)[[1]]

Process raw count data

Ferr <- Ferr[rowSums(is.na(Ferr)) == 0, ]

countdata_Fer <- countdata[ferr_genes_4plots,]  # select genes here

Fer_samples = colsplit(colnames(countdata_Fer), pattern = "_", names = c("Population", "Time", "Replicate"))
Fer_plot = cbind(Fer_samples, t(countdata_Fer))
Fer_melt = melt(data = Fer_plot, id.vars = c("Population", "Time", "Replicate"), measure.vars = ferr_genes_4plots)  # select genes here
Fer_dat = Rmisc::summarySE(Fer_melt, measurevar = "value", groupvars = c("Population", "Time", "variable"))
Fer_dat$Time = as.numeric(gsub("[^[:digit:]]","",Fer_dat$Time))

Plot change over time of BRAFi for each subclone

Changed plotting from linear to log2 scale for better visualization of changes.

Fer_ggploted <- ggplot(Fer_dat, aes(x=Time, y=log2(value), group = interaction(variable, Population))) + 
  geom_line(linewidth=1.5, aes(color = Population)) + 
  geom_point(size = 1.5, aes(color = Population)) + facet_wrap(~variable, ncol = 5, scales = "free") +
  geom_errorbar(aes(ymin=log2(value-sd), ymax=log2(value+sd), color = Population), width=.2, linewidth=1.5) +
  theme_bw() + xlab("Time (days)") + ylab("Gene Counts") +
  ggtitle("Ferroptosis gene signature") +
  theme(legend.text = element_text(size = 10), legend.position = "right", 
        plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), axis.text=element_text(size=12),
        legend.title = element_text(size=12,face="bold"),
        axis.title=element_text(size=12,face="bold"))

Fer_ggploted

if(SAVEPLOTS) ggsave("FerroptosisGeneSignature_rawCounts_SKMEL5sublines+treatment.pdf", width = 20, height = 25)

NOTE: This will overwrite objects with new data

Normalized data from DESeq2

normdata_Fer <- as.data.frame(assay(rld))[ferr_genes_4plots,]

Fer_match <- colsplit(colnames(normdata_Fer), pattern = "_", names = c("Population", "Time", "Replicate"))
Fer_plot <- cbind(Fer_match, t(normdata_Fer))
Fer_melt <- melt(data = Fer_plot, id.vars = c("Population", "Time", "Replicate"), measure.vars = unique(colnames(Fer_plot))[4:42])
Fer_dat <- Rmisc::summarySE(Fer_melt, measurevar = "value", groupvars = c("Population", "Time", "variable"))
Fer_dat$Time <- as.numeric(gsub("[^[:digit:]]","",Fer_dat$Time))

# normdata_Fer
samples <- c("SC01_day0", "SC01_day3", "SC01_day8", "SC07_day0", "SC07_day3", "SC07_day8", "SC10_day0", "SC10_day3", "SC10_day8")
test <- normdata_Fer
test1 <- sapply(samples, function(x) rowMeans(test[, grep(x, colnames(test))]))
test2 <- as.data.frame(test1[complete.cases(test1),])

####
# Plot pheatmap (z-score) of gene list 
test3 <- test2
test3$Gene <- rownames(test2)
test3_sub <- subset(test3, !Gene %in% excl_genes)
test3_sub <- test3_sub[,1:9]
pheatmap(test3_sub,cluster_cols=FALSE, cluster_rows = F, scale = "row")

test4_special <- subset(test3, Gene %in% ferr_genes_4plots)
test4_special <- test4_special[,1:9]
pheatmap(test4_special,cluster_cols=FALSE, cluster_rows = F, scale = "row")

NA
NA

allFC <- function(DEProc,startcol,endcol){ 
  GE_fold = DEProc[,-c(startcol:endcol)]
  colvec = colnames(DEProc)[startcol:endcol]
  
  #Last index is a self comparison and is removed
  for(k in 1:(length(colvec)-1)){
    #Start with column that is 1 away from index 
    for(j in (k+1):length(colvec)){
      compnam = paste0(colvec[j],"/",colvec[k])
      #Loop through each gene/row  
      for(i in 1:nrow(DEProc)){
        f = DEProc[i,colvec[j]]
        h = DEProc[i,colvec[k]]
        GE_fold[i, compnam] = log2(f/h)

        #Capture upregulation and down regulation
        # if(f>h){
        #   # GE_fold[i,compnam] = 2^(f-h)
        #   GE_fold[i,compnam] = log2(f/h)
        # }else{
        #   # GE_fold[i,compnam] = -2^(h-f)
        #   GE_fold[i,compnam] = log2(f/h)
      #   }
      #   
     }
    }
  }
  
  return(GE_fold)
  
}

GE_fold <- allFC(test2, 1,9)
ImpRat = c("SC01_day3/SC01_day0", "SC01_day8/SC01_day0", 
           "SC07_day3/SC07_day0", "SC07_day8/SC07_day0", 
           "SC10_day3/SC10_day0", "SC10_day8/SC10_day0")
Imp_fold <- GE_fold[,ImpRat]

Ferr <- read.delim("../data/KEGGFerroptosis_hsa04216_06-25-18.txt", header=T, stringsAsFactors = F)
Ferr_fold <- Imp_fold[rownames(Imp_fold) %in% ferr_genes_4plots,]

Ferr_fold$Gene <- rownames(Ferr_fold)

# clean column names
colnames(Ferr_fold) <- gsub("/SC[01][170]_day0","",colnames(Ferr_fold))

Ferr_fold_melt <- melt(Ferr_fold, id.vars = "Gene")
Ferr_fold_melt$Gene <- factor(Ferr_fold_melt$Gene, levels = rev(ferr_genes_4plots))

Ferr_fold_melt <- subset(Ferr_fold_melt, !Gene %in% excl_genes)

myplot <- ggplot(Ferr_fold_melt, aes(variable, Gene, fill = value)) + 
  geom_tile(color = "black") + theme_bw() +
  scale_fill_gradientn(
    colors=c("blue","white","red","red4"),
    values=rescale(c(min(Ferr_fold_melt$value), 0, max(Ferr_fold_melt$value)/2, max(Ferr_fold_melt$value))),
    limits=c(min(Ferr_fold_melt$value),max(Ferr_fold_melt$value)),
    name = "Log2 Fold Change"
  ) + ylab("") + xlab("") + 
  theme(axis.text=element_text(size=10),
        axis.text.x=element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=12),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

myplot

if(SAVEPLOTS) ggsave("Ferroptosis_FCHM_selected_wide.pdf", width = 8, height = 8)

Using ComplexHeatmap()

if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("ComplexHeatmap")